package com.asa.numanaly;

import com.asa.HanShu;

/**
 * 常微分方程的边值问题 我很期待这个和I里面的逼近特征值一样 这个可以用在刚体里面
 * 第五章也就是E里面也有求常微分方程初值问题，但那个只能是一阶的，而且只满足一阶条件
 * 
 * @author asa
 *
 */
public class K {
	
	
	public static void main(String[] args) {
		double a = 1;
		double b = 2;
		double aa = 1;
		double bb = 2;
		int N = 9;

		HanShu hanshu2 = new HanShu() {

			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double x0) {
				// TODO Auto-generated method stub
				return -2.0/x0;
			}
		};
		HanShu hanshu1 = new HanShu() {

			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double x0) {
				// TODO Auto-generated method stub
				return 2.0/(x0*x0);
			}
		};
		HanShu hanshu0 = new HanShu() {

			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double x0) {
				// TODO Auto-generated method stub
				return Math.sin(Math.log(x0))/(x0*x0);
			}
		};

		changweifenfangcheng3(a, b, aa, bb, N, hanshu2, hanshu1, hanshu0);
	}

	/*
	 * p591 先实现这个，这个我读大学的时候实现过
	 * 
	 */
	private static void changweifenfangcheng3(double a, double b, double aa, double bb, int N, HanShu hanshu2,
			HanShu hanshu1, HanShu hanshu0) {
		double h = (b-a)/(N+1);
		double x = a+h;
		
		double[] ai =new double[N+1];
		double[] bi =new double[N+1];
		double[] ci =new double[N+1];

		double[] di =new double[N+1];

		
		ai[1] = 2+h*h*hanshu1.hanshu(x);
		bi[1] = -1+(h/2)*hanshu2.hanshu(x);
		di[1] = -h*h*hanshu0.hanshu(x)+ (1+(h/2)*hanshu2.hanshu(x))*aa;
		
		
		for (int i = 2; i < N; i++) {
			
			x = a+i*h;
			ai[i] = 2+h*h*hanshu1.hanshu(x);
			bi[i] = -1 + (h/2)*hanshu2.hanshu(x);
			ci[i] = -1 - (h/2)*hanshu2.hanshu(x);
			di[i] = -h*h*hanshu0.hanshu(x);
			
		}
		
		
		x = b-h;
		ai[N] = 2+h*h*hanshu1.hanshu(x);
		ci[N] = -1-(h/2)*hanshu2.hanshu(x);
		di[N] = -h*h*hanshu0.hanshu(x)+(1-(h/2)*hanshu2.hanshu(x))*bb;
		
		
		
		
		
		
		double[] l =new double[N+1];
		double[] u =new double[N+1];
		double[] z =new double[N+1];
		l[1] = ai[1];
		u[1] = bi[1]/ai[1];
		z[1] = di[1]/l[1];
		
		for (int i = 2; i <= N-1; i++) {
			l[i] = ai[i] - ci[i]*u[i-1];
			u[i] = bi[i]/l[i];
			z[i] = (di[i]-ci[i]*z[i-1])/l[i];
			
			
		}
		l[N] = ai[N]-ci[N]*u[N-1];
		z[N] = (di[N]-ci[N]*z[N-1])/l[N];
		
		double[] w = new double[N+2];
		w[0] = aa;
		w[N+1] = bb;
		w[N] = z[N];
		for (int i = N-1; i >0; i--) {
			w[i] = z[i]-u[i]*w[i+1];
		}
		for (int i = 0; i <= N+1; i++) {
			x = a+i*h;
			System.out.println(x+"          "+w[i]);
		}
	}
	

	

}
